/* This file has the declaration for Riemann*Tensor4 yielding
   a double.  I simplify the expression by removing the identically
   zero components of Riemann. */

#pragma once

namespace FTensor
{
  /* A(i,j,k,l)*B(i,j,k,l) */

  template <class A, class B, char i, char j, char k, char l>
  double operator*(const Riemann_Expr<A, i, j, k, l> &a,
                   const Tensor4_Expr<B, i, j, k, l> &b)
  {
    return a(1, 0, 0, 1) * b(1, 0, 0, 1) + a(2, 0, 0, 1) * b(2, 0, 0, 1)
           + a(0, 1, 0, 1) * b(0, 1, 0, 1) + a(2, 1, 0, 1) * b(2, 1, 0, 1)
           + a(0, 2, 0, 1) * b(0, 2, 0, 1) + a(1, 2, 0, 1) * b(1, 2, 0, 1)
           + a(1, 0, 0, 2) * b(1, 0, 0, 2) + a(2, 0, 0, 2) * b(2, 0, 0, 2)
           + a(0, 1, 0, 2) * b(0, 1, 0, 2) + a(2, 1, 0, 2) * b(2, 1, 0, 2)
           + a(0, 2, 0, 2) * b(0, 2, 0, 2) + a(1, 2, 0, 2) * b(1, 2, 0, 2)
           + a(1, 0, 1, 0) * b(1, 0, 1, 0) + a(2, 0, 1, 0) * b(2, 0, 1, 0)
           + a(0, 1, 1, 0) * b(0, 1, 1, 0) + a(2, 1, 1, 0) * b(2, 1, 1, 0)
           + a(0, 2, 1, 0) * b(0, 2, 1, 0) + a(1, 2, 1, 0) * b(1, 2, 1, 0)
           + a(1, 0, 1, 2) * b(1, 0, 1, 2) + a(2, 0, 1, 2) * b(2, 0, 1, 2)
           + a(0, 1, 1, 2) * b(0, 1, 1, 2) + a(2, 1, 1, 2) * b(2, 1, 1, 2)
           + a(0, 2, 1, 2) * b(0, 2, 1, 2) + a(1, 2, 1, 2) * b(1, 2, 1, 2)
           + a(1, 0, 2, 0) * b(1, 0, 2, 0) + a(2, 0, 2, 0) * b(2, 0, 2, 0)
           + a(0, 1, 2, 0) * b(0, 1, 2, 0) + a(2, 1, 2, 0) * b(2, 1, 2, 0)
           + a(0, 2, 2, 0) * b(0, 2, 2, 0) + a(1, 2, 2, 0) * b(1, 2, 2, 0)
           + a(1, 0, 2, 1) * b(1, 0, 2, 1) + a(2, 0, 2, 1) * b(2, 0, 2, 1)
           + a(0, 1, 2, 1) * b(0, 1, 2, 1) + a(2, 1, 2, 1) * b(2, 1, 2, 1)
           + a(0, 2, 2, 1) * b(0, 2, 2, 1) + a(1, 2, 2, 1) * b(1, 2, 2, 1);

    //    return a(0,0,0,0)*b(0,0,0,0) + a(1,0,0,0)*b(1,0,0,0) +
    //    a(2,0,0,0)*b(2,0,0,0)
    //      + a(0,1,0,0)*b(0,1,0,0) + a(1,1,0,0)*b(1,1,0,0) +
    //      a(2,1,0,0)*b(2,1,0,0)
    //      + a(0,2,0,0)*b(0,2,0,0) + a(1,2,0,0)*b(1,2,0,0) +
    //      a(2,2,0,0)*b(2,2,0,0)
    //      + a(0,0,0,1)*b(0,0,0,1) + a(1,0,0,1)*b(1,0,0,1) +
    //      a(2,0,0,1)*b(2,0,0,1)
    //      + a(0,1,0,1)*b(0,1,0,1) + a(1,1,0,1)*b(1,1,0,1) +
    //      a(2,1,0,1)*b(2,1,0,1)
    //      + a(0,2,0,1)*b(0,2,0,1) + a(1,2,0,1)*b(1,2,0,1) +
    //      a(2,2,0,1)*b(2,2,0,1)
    //      + a(0,0,0,2)*b(0,0,0,2) + a(1,0,0,2)*b(1,0,0,2) +
    //      a(2,0,0,2)*b(2,0,0,2)
    //      + a(0,1,0,2)*b(0,1,0,2) + a(1,1,0,2)*b(1,1,0,2) +
    //      a(2,1,0,2)*b(2,1,0,2)
    //      + a(0,2,0,2)*b(0,2,0,2) + a(1,2,0,2)*b(1,2,0,2) +
    //      a(2,2,0,2)*b(2,2,0,2)
    //      + a(0,0,1,0)*b(0,0,1,0) + a(1,0,1,0)*b(1,0,1,0) +
    //      a(2,0,1,0)*b(2,0,1,0)
    //      + a(0,1,1,0)*b(0,1,1,0) + a(1,1,1,0)*b(1,1,1,0) +
    //      a(2,1,1,0)*b(2,1,1,0)
    //      + a(0,2,1,0)*b(0,2,1,0) + a(1,2,1,0)*b(1,2,1,0) +
    //      a(2,2,1,0)*b(2,2,1,0)
    //      + a(0,0,1,1)*b(0,0,1,1) + a(1,0,1,1)*b(1,0,1,1) +
    //      a(2,0,1,1)*b(2,0,1,1)
    //      + a(0,1,1,1)*b(0,1,1,1) + a(1,1,1,1)*b(1,1,1,1) +
    //      a(2,1,1,1)*b(2,1,1,1)
    //      + a(0,2,1,1)*b(0,2,1,1) + a(1,2,1,1)*b(1,2,1,1) +
    //      a(2,2,1,1)*b(2,2,1,1)
    //      + a(0,0,1,2)*b(0,0,1,2) + a(1,0,1,2)*b(1,0,1,2) +
    //      a(2,0,1,2)*b(2,0,1,2)
    //      + a(0,1,1,2)*b(0,1,1,2) + a(1,1,1,2)*b(1,1,1,2) +
    //      a(2,1,1,2)*b(2,1,1,2)
    //      + a(0,2,1,2)*b(0,2,1,2) + a(1,2,1,2)*b(1,2,1,2) +
    //      a(2,2,1,2)*b(2,2,1,2)
    //      + a(0,0,2,0)*b(0,0,2,0) + a(1,0,2,0)*b(1,0,2,0) +
    //      a(2,0,2,0)*b(2,0,2,0)
    //      + a(0,1,2,0)*b(0,1,2,0) + a(1,1,2,0)*b(1,1,2,0) +
    //      a(2,1,2,0)*b(2,1,2,0)
    //      + a(0,2,2,0)*b(0,2,2,0) + a(1,2,2,0)*b(1,2,2,0) +
    //      a(2,2,2,0)*b(2,2,2,0)
    //      + a(0,0,2,1)*b(0,0,2,1) + a(1,0,2,1)*b(1,0,2,1) +
    //      a(2,0,2,1)*b(2,0,2,1)
    //      + a(0,1,2,1)*b(0,1,2,1) + a(1,1,2,1)*b(1,1,2,1) +
    //      a(2,1,2,1)*b(2,1,2,1)
    //      + a(0,2,2,1)*b(0,2,2,1) + a(1,2,2,1)*b(1,2,2,1) +
    //      a(2,2,2,1)*b(2,2,2,1)
    //      + a(0,0,2,2)*b(0,0,2,2) + a(1,0,2,2)*b(1,0,2,2) +
    //      a(2,0,2,2)*b(2,0,2,2)
    //      + a(0,1,2,2)*b(0,1,2,2) + a(1,1,2,2)*b(1,1,2,2) +
    //      a(2,1,2,2)*b(2,1,2,2)
    //      + a(0,2,2,2)*b(0,2,2,2) + a(1,2,2,2)*b(1,2,2,2) +
    //      a(2,2,2,2)*b(2,2,2,2);
  }

  template <class A, class B, char i, char j, char k, char l>
  double operator*(const Tensor4_Expr<B, i, j, k, l> &b,
                   const Riemann_Expr<A, i, j, k, l> &a)

  {
    return operator*(a, b);
  }
}
